Variational RVB wave function for the spin- 1/2 Heisenberg Model on honeycomb lattice 
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In this work, a long-range resonating valence bond state is proposed as a variational wave function for the 
ground state of the S — 1/2 antiferromagnetic Heisenberg model on the honeycomb lattice. Employing Vari- 
ational Monte Carlo (VMC) method, we show that the ground state energy obtained from such RVB wave 
function, lies well below the energy of the Neel state and compares very well to the energies evaluated from 
spin-wave theory and series expansion method. We also obtain the spin-spin correlation function along zig-zag 
and armchair directions and find that the two correlations are different, which indicates the anisotropic nature of 
the system. We compare our results with the square lattice and we show that although the quantum fluctuations 
on honeycomb lattice are much stronger, but do completely not destroy the magnetic order. 

PACS numbers: 75.10.Jm, 75.50.Ee 



Introduction: The resonating valence bond (RVB) state was 
originally proposed by Anderson and Fazekas as a possible 
ground state for the S = 1/2 Heisenberg spins with near- 
est neighbor anti-ferromagnetic coupling on triangular lat- 
tice iQllII]- They found that for anisotropic Heisenberg model, 
near the Ising limit, the liquid-like RVB state is energeti- 
cally more favorable than the Neel state. In such a state 
the S = 1/2 atoms residing on the lattice points, form sin- 
glet valence bonds in pairs and so lose some of their anti- 
ferromagnetic exchange energy with respect to the Neel or- 
der. In order to regain some of the lost exchange energy, they 
have to resonate quantum mechanically among many different 
pairing configurations. These states form the basis of Paul- 
ing's early theories of aromatic molecules such as benzene fl, 
however his theory was unable to give a proper description of 
the metallic state. The idea that the RVB state of spin-liquid 
type may give a precise picture of the two dimensional quan- 
tum anti-ferromagnet was once the most attractive topic after 
the discovery of the high-T c superconductivity, when Ander- 
son |01 suggested that an RVB state naturally leads to incipient 
superconductivity from preformed singlet pairs in the parent 
insulating state. 

The anti-ferromagnetic Heisenberg model arises naturally 
in a Mott insulator, in which a system with an odd number of 
electrons per unit cell is insulating due to the strong Coulomb 
repulsion between two electrons on the same site ([/). In such 
systems the kinetic exchange mechanism due to virtual hop- 
ping between anti-parallel spin configurations leads to anti- 
ferromagntic exchange J = 4t 2 /U between the spins, where 
t is the hopping integral ]5|]. For S = 1/2 and in low di- 
mensions, strong quantum fluctuations make the RVB liquid 
more favorable than the classically ordered Neel state. RVB 
is a fair description of anti-ferromagnetic Heisenberg linear 
chain, where there is no long-range order according to the 



exact Bethe's solution The S = 1/2 anti-ferromagnet 
Heisenberg model on square lattice has been extensively stud- 
ied |0] and it has been established that there is a Neel order 
in the system, although the quantum fluctuations reduce the 
staggered magnetization with respect to its classical value. 

Apart from the fabrication of Graphene sheets which 
brought the honeycomb structure to the focus of condensed 
matter community, the recent discovery of compounds such 
as InCu 2 /3Vi/303 H in which the Cu +2 ions in the 
copper-oxide layers form a two-dimensional S = 1/2 anti- 
ferromagnt Heisenberg on a honeycomb lattice are our mo- 
tivations for this study. Also the recent progress in the field 
of ultracold atoms and trapping techniques |9[] along with the 
ability to tune the interaction parameters via the Feshbach res- 
onance ifioll can be thought of another way to realize Heisen- 
berg spins (of localized fermions) on a honeycomb optical lat- 
tice. 

In this work, we study the 5=1/2 Heisenberg spins with 
nearest-neighbor anti-ferromagnetic interactions on the hon- 
eycomb lattice. The number of nearest neighbors in honey- 
comb lattice equals 3, which is less than 4 of the square lat- 
tice, leading to enhancement of quantum fluctuations. This 
suggests that the RVB state could be a better variational 
wave function for honeycomb than the square lattice. In 
this work we choose the variational RVB wave function pro- 
posed by Liang, Doucot and Anderson in the context of HTSC 
cuprates 1 1 ill and show that RVB ansatz in honeycomb lattice 
gives very good results for the ground state energy by com- 
paring with other methods. 

Numerical calculation: Consider an system of atoms with 
S = l/2onaa honeycomb lattice, consisting of two interpen- 
etrating bravais sublattices (Fig. [TJ. The Heisenberg Hamilto- 
nian with anti-ferromagnetic nearest neighbor interactions for 
this system is given by: 
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where J > and (i, j) denotes the nearest neighbors. Since 
the ground state of such Hamiltonian is a spin singlet 11211 . we 
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FIG. 1: Honeycomb lattice as a superposition of two simple mono- 
clinic Bravais lattices denoted by filled and empty circles. 



employ a RVB state as trial wave function. After its optimiza- 
tion, we calculate the ground state properties. The RVB trial 
wave function can be considered as summation over the all 
possible singlet pairings between each spin in one sublattice 
(say A) to spins on the other sublattice (B). Therefore the trial 
wave function is given by: 

|*> = h ^ il -h)- h (in ~ jn){il,jl)-(in,3n) 

= yi"Kct)ict) (2) 

i 

Where represent a singlet bond. Denoting spin up state 
(| T)) by a and spin down state (| |)) by (3, it can be expressed 
as: 

= (Oj/?j - PiOtj). 



(i,3) 



(3) 



In Eq. (O, \c) stands for a valence bond configuration which 
can be represented by: 



\c) = U(la,ja). (4) 

£1=1 

The weight of a configuration |c) of valence bonds is given by 

n 

w(c) = Y[ h(i a - ja), (5) 

a=l 

where h is a pairing function describing a singlet bond as a 
function of bond length. 

This wave function contains the two limiting cases of the 
nearest-neighbor RVB liquid (when h(l) = for I > 1) 
and the Neel state (for which h(l) is independent of the bond 
length). Sutherland derived a set of simple rules for esti- 
mating (ci|Sj.Sj|c2) 01311 . It is easy to show that the over- 
lap between two configurations |ci) and \c 2 ) is given by 




FIG. 2: The loop covering (c\\c2) is the superposition of |ci) (solid 
line) and \c2) (dashed line) on a honeycomb lattice. |ci) and |c2) are 
two singlet valence bond configurations with equal weight factors. 
In this example, there are six loops. 



(ci|c2) = 2 N ( Cl ' C2 \ where N(c\, c 2 ) is the number of loops 
in the overlap of the two configurations (Fig. |2j. To compute 
the ground state energy as well as the spin-spin correlation 
functions, we use the following rules for the matrix elements 
of the Hamiltonian: (i) (ci|S.;.Sj \c 2 ) = 0, if i and j belong to 



two different loops; (ii) (ci|Si.Sj|c2) = ±|(ci|c2), if i and 
j belong to the same loop, with a minus sign when i and j 
belong to two different sublattices, and a plus sign otherwise. 

Using the wave function Q along with the Hamiltonian (Q~|), 
the ground state energy is given by: 



mm 

E w(ci)w(c 2 )(ci\c 2 ) (ci\H\c 2 ) 
(*|*> X ( C1 |C 2 > 

Ci,C2 vi/ \i/ 

^P(ci,c 2 ) x E( Cl ,c 2 ), (6) 



in which E(c\ ,c 2 ) is the contribution in the ground states en- 
ergy arising from the overlap of configurations \c\) and \c 2 ) 
with the weight P(c lf c 2 ) = w{ci)w { ^ { } cilc2) . Such overlaps 
can be graphically represented by loop coverings depicted in 
Fig. [2l 1140 . The number of terms in this expression diverges 
exponentially with the system size, making the exact evalua- 
tion of the summation impossible. To overcome this difficulty, 
we use Monte Carlo approach based on important sampling 
for evaluating the ground state energy and also the spin-spin 
correlations, since P(ci, c 2 ) is positive, here there is no mi- 
nus sign problem, and the Monte Carlo estimates can be made 
very precise. 

Consider two configurations c\ and c 2 with a given loop 
coverage. Using the standard Metropolis algorithm, loop con- 
figurations are updated by randomly choosing a couple of sites 
and exchanging their end point connections with a probability 
that satisfies the detailed balance condition 111 111 . The matrix 
elements are evaluated according to Sutherland's rule, once 
the loop covering associated with the two configurations is 
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FIG. 3: The spin-spin correlation as a function of distance for the 
state labeled (1, 3, 5). See the text for explanations of the state. Cir- 
cles denote zigzag direction data and triangles denote armchair data. 



known. In our calculations, boundary conditions are take to 
be periodic. The equilibrium state does not depend on the 
initial state, but in order to reach the equilibrium distribution 
faster, we started with a dimer state. In following we discuss 
our results for short-range and long-range RVB states. 

Short-range RVB wave function: First we consider some 
wave functions with short singlet bonds. One example of this 
type is the nearest-neighbor RVB (NNRVB) trial wave func- 
tions for which h(l) = 1 and h(l) = for larger distances 
I > 1. We choose I to be the Manhatan distance defined 
for two points separated by a vector n u u + n v v in Fig. [T|by 
I — \n u \ + \n v \. Defining a; = j^'-i] ' we investigate other 



h(2l- 

examples such as (1,3) state with one variational parameter 
a\ = and a/ = for I > 2; (1,3,5) state with two 

variational parameters ai = 02 = ^§y; and exponential 
state for which h{l) decays exponentially for large distances. 
In this state we have three variational parameters 01, a 2 and 
ai = const, for I > 2. Among the short range states, the ex- 
ponential state has lowest energy (see TableQ]i. The statistical 
error in evaluation of energy from Monte Carlo calculation is 
0.0002 and we check that the finite-size effects are less than 
this error. We have also computed the spin-spin correlation 
function defined as: 



Cij = ^P(ci,c 2 ) 

Ci,C2 



(ci|Sj.Sj|c 2 ) 

<C1|C2> 



(7) 



which results in exponentially decaying of the correlation at 
large distances, indicating the absence of long- range order for 
this variational states (Fig. [3). As it can be seen in Fig. [3] the 
correlation length along zigzag direction is larger than arm- 
chair direction. The correlation lengths £ in units of the lattice 
spacing are listed in Table U 

Long-range wave functions: Next we study spin-spin cor- 
relations and energy for the class of wave functions which be- 
have like h(l) oc l~ p in long distances. To optimize the en- 



Hi) 
h(i) = 1 

h(l) =0, / > 1 

01 = § 
ai = 0, I > 1 



NNRVB 0.4310(2) £ a = 0.84 
£r = l 

(1, 3) 0.5087(2) £ a = 4.7 
£z = 5.9 

(1,3,5) 0.5319(2) £ a = 9.55 (01, a 2 ) = (i, f) 
£ 2 = 11.1 h(l) = 0, / > 2 
Exponential 0.5437(2) (a 1; a 2 ) = (£, |) 
ai = 0.32, I > 3 

TABLE I: Ground state energy, spin-spin correlation length and op- 
timized parameters for the short ranged RVB states. The subscripts 
a, z stand for armchair and zigzag directs, respectively. 
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FIG. 4: spin-spin correlation as a function of distance for the states 
with algebraically decaying states with p = 3 (see the text), versus r 
along the zigzag (circle) and armchair (triangle) directions. 



ergy, we choose a\ 



h (3) 
h{V> 



, a 3 



h(5) 
h(3) 



as variational param- 



eters and optimize them for various exponents p defined by 
h(l) = h(5)(5/l) p for I > 5. For long-bond singlets the spin- 
spin correlation function does not vanish at large distances and 
so the staggered magnetization can be determined from the 
tails of the correlation function. Because of the anisotropic 
nature of the honeycomb lattice, the correlations along zigzag 
and armchair directions are different hence we define the mag- 
netization by 



Czigzagfr) "l" C armc h a i r (r) 



M s = J lim 

Y 2 

where L is the linear size of the system. 



(8) 



Ms (ai,a 2 ) 



P 

2 0.5434(2) 0.30(2) 
2.5 0.5440(2) 0.26(2) 

3 0.5438(2) 0.25(2) 
3.5 0.5431(2) 0.20(2) 

4 0.5430(2) 0.19(2) 



(— -) 
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TABLE II: Ground state energy, staggered magnetization and opti- 
mized parameters for the long ranged RVB states. 
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The optimal values for the parameters as well as ground 
state energies and staggered magnetizations within this class 
of states are listed in Table [II] for various values p = 
2, 2.5, 3, 3.5 and 4. For this set of exponents, we found that 
energy is minimum for p = 2.5 with E a = —0.5440 J ± 
0.0002J, which is lower than the energy of the short 
ranged RVB states. The staggered magnetization for p = 
2, 2.5, 3, 3.5 and 4 are %60, %52, %50, %40, and %38 of Neel 
state, respectively. The magnetization decreases when p in- 
creases and seems to disappear when p > 4. 

Conclusion: In summary, we presented a variational Monte 
Carlo estimate of the ground state energy for the S = 1/2 
AF Heisenberg model on honeycomb lattice employing both 
short-range and l ong -range RVB trial wave functions. Similar 
to square lattice ll ill , we found that the long -bond wave func- 
tions give a better description of the ground state in honey- 
comb lattice. The optimized values of ground state energy and 
staggered magnetization are E /J = —0.5440 ± 0.0002 and 
M s = 0.26, respectively. The ratio of the RVB ground state 
energy to the Neel state energy is 



0.5540 



1., 



E N 0.375 

RVB energy is %48 less than the Neel state. In square lattice, 
this ratio is %34, which can be justified in terms of stronger 
quantum fluctuation in honeycomb lattice due to its smaller 
coordination number. The quantum fluctuations also reduce 
the order parameter by about %50 with respect to classical 
Neel order. 

Our results are in excellent agreement with other numerical 
methods, such as the series expansion lfl~5ll which gives the 
ground state energy ^f- = —0.5443 ± 0.0003 and magnetiza- 
tion M s = 0.266 ± 0.009 as well as spin wave calculation in 



first approximation with 



-0.5324 and M s 



0.2418 



and in second approximation with ^j- = —0.5489 and M s = 
0.2418 0. This shows that RVB picture captures the main 
properties of the ground state for the AF Heisenberg model 
on the honeycomb lattice with higher precision than square 
lattice, leading to the conclusion that the stronger are quan- 
tum fluctuations, the more favorable is the RVB ground state. 
Therefore it seems that the RVB state is a good candidate for 
the ground state of the AF Heisenberg model on the honey- 
comb lattice. One can take a long ranged RVB as a reference 
starting point to proceed with the calculations of excitation 
energies lfl7ll . Including charge fluctuations in such an RVB 



state by addi ng hopp ing of underlying electrons has also been 
investigated I18U19I1 . This view point can be taken as a possi- 
ble rout to describe the Dirac liquid semi-metallic in terms of 
RVB wave functions 1201. 
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